if (!require("interp")) {install.packages("interp")}
library(ISLR2) # Para el conjunto de datos Wage
library(splines) # Para splines naturales (ns)
library(gam) # Para la función gam()14 Modelos Generalizados Aditivos
15 Introducción
Este documento explora los Modelos Aditivos Generalizados (GAM), una extensión flexible de los modelos lineales que permite modelar la relación entre las variables predictoras y la respuesta mediante funciones de suavizado no lineales. A diferencia de un modelo lineal estándar que asume una relación lineal, un GAM modela la respuesta como una suma de funciones suaves de los predictores.
El modelo tiene la forma: \[ E(Y|X_1, ..., X_p) = \beta_0 + f_1(X_1) + f_2(X_2) + \dots + f_p(X_p) \] Donde \(f_j()\) son funciones de suavizado no especificadas (como splines o regresiones locales) que se estiman a partir de los datos.
15.0.1 Librerías Requeridas
15.0.2 Carga de Datos
Se utiliza el conjunto de datos Wage de la librería ISLR2.
data(Wage)
attach(Wage) # Permite llamar a las columnas por su nombre directamente
# View(Wage)16 a) Modelos Aditivos Generalizados (GAM)
16.1 Estimación de un GAM con lm() y Splines Naturales
Un GAM puede ser aproximado usando un modelo lineal (lm()) donde las variables no lineales se modelan mediante splines naturales (ns()). Los splines son funciones polinómicas por tramos que se unen suavemente en puntos llamados “nodos”.
# Se modela 'year' y 'age' de forma no paramétrica con splines naturales,
# y 'education' de forma paramétrica.
gam1 <- lm(wage ~ ns(year, df = 4) + ns(age, df = 5) + education, data = Wage)
summary(gam1)
Call:
lm(formula = wage ~ ns(year, df = 4) + ns(age, df = 5) + education,
data = Wage)
Residuals:
Min 1Q Median 3Q Max
-120.513 -19.608 -3.583 14.112 214.535
Coefficients:
Estimate Std. Error t value Pr(>|t|)
(Intercept) 46.949 4.704 9.980 < 2e-16 ***
ns(year, df = 4)1 8.625 3.466 2.488 0.01289 *
ns(year, df = 4)2 3.762 2.959 1.271 0.20369
ns(year, df = 4)3 8.127 4.211 1.930 0.05375 .
ns(year, df = 4)4 6.806 2.397 2.840 0.00455 **
ns(age, df = 5)1 45.170 4.193 10.771 < 2e-16 ***
ns(age, df = 5)2 38.450 5.076 7.575 4.78e-14 ***
ns(age, df = 5)3 34.239 4.383 7.813 7.69e-15 ***
ns(age, df = 5)4 48.678 10.572 4.605 4.31e-06 ***
ns(age, df = 5)5 6.557 8.367 0.784 0.43328
education2. HS Grad 10.983 2.430 4.520 6.43e-06 ***
education3. Some College 23.473 2.562 9.163 < 2e-16 ***
education4. College Grad 38.314 2.547 15.042 < 2e-16 ***
education5. Advanced Degree 62.554 2.761 22.654 < 2e-16 ***
---
Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
Residual standard error: 35.16 on 2986 degrees of freedom
Multiple R-squared: 0.293, Adjusted R-squared: 0.2899
F-statistic: 95.2 on 13 and 2986 DF, p-value: < 2.2e-16Nota: En la salida de summary(gam1), solo los coeficientes de las variables paramétricas (como education) son directamente interpretables. Los coeficientes asociados a los splines (ns(...)) no tienen una interpretación sencilla.
16.2 Estimación con la Librería gam
La forma recomendada de ajustar un GAM es con la librería gam, que está diseñada específicamente para este propósito. Utiliza la función s() para especificar términos de suavizado (en este caso, smoothing splines).
# s(variable, df) define un término de suavizado con 'df' grados de libertad.
gam.m3 <- gam(wage ~ s(year, 4) + s(age, 5) + education, data = Wage)
summary(gam.m3)
Call: gam(formula = wage ~ s(year, 4) + s(age, 5) + education, data = Wage)
Deviance Residuals:
Min 1Q Median 3Q Max
-119.43 -19.70 -3.33 14.17 213.48
(Dispersion Parameter for gaussian family taken to be 1235.69)
Null Deviance: 5222086 on 2999 degrees of freedom
Residual Deviance: 3689770 on 2986 degrees of freedom
AIC: 29887.75
Number of Local Scoring Iterations: NA
Anova for Parametric Effects
Df Sum Sq Mean Sq F value Pr(>F)
s(year, 4) 1 27162 27162 21.981 2.877e-06 ***
s(age, 5) 1 195338 195338 158.081 < 2.2e-16 ***
education 4 1069726 267432 216.423 < 2.2e-16 ***
Residuals 2986 3689770 1236
---
Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
Anova for Nonparametric Effects
Npar Df Npar F Pr(F)
(Intercept)
s(year, 4) 3 1.086 0.3537
s(age, 5) 4 32.380 <2e-16 ***
education
---
Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1Nota: La salida de summary.gam no muestra los coeficientes de los términos no paramétricos, ya que no son el foco de la interpretación. En su lugar, proporciona una tabla de ANOVA para evaluar la significancia de cada término.
17 b) Interpretación Gráfica del GAM
La principal herramienta para interpretar un GAM es la visualización de las funciones de suavizado estimadas para cada predictor. Los gráficos muestran la forma de la relación entre cada predictor y la respuesta, manteniendo constantes las demás variables.
par(mfrow = c(1, 3))
# Gráfico de los efectos parciales del modelo ajustado con gam()
plot(gam.m3, se = TRUE, col = "blue")
18 c) Comparación de Modelos con anova()
Se utiliza la prueba F para comparar formalmente la bondad de ajuste de modelos anidados con diferentes especificaciones para la variable year.
- Modelo 1: Sin
year. - Modelo 2: Con
yearcomo término lineal (paramétrico). - Modelo 3: Con
yearcomo término de suavizado (no paramétrico).
gam.m1 <- gam(wage ~ s(age, 5) + education, data = Wage)
gam.m2 <- gam(wage ~ year + s(age, 5) + education, data = Wage)
# Prueba F para comparar los tres modelos anidados
anova(gam.m1, gam.m2, gam.m3, test = "F")Analysis of Deviance Table
Model 1: wage ~ s(age, 5) + education
Model 2: wage ~ year + s(age, 5) + education
Model 3: wage ~ s(year, 4) + s(age, 5) + education
Resid. Df Resid. Dev Df Deviance F Pr(>F)
1 2990 3711731
2 2989 3693842 1 17889.2 14.4771 0.0001447 ***
3 2986 3689770 3 4071.1 1.0982 0.3485661
---
Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1Interpretación del anova(): - La primera fila (p-valor < 0.001) compara el Modelo 1 con el Modelo 2. Se rechaza \(H_0\), lo que significa que incluir year de forma lineal mejora significativamente el modelo. - La segunda fila (p-valor = 0.348) compara el Modelo 2 con el Modelo 3. No se rechaza \(H_0\), lo que sugiere que modelar year con una función de suavizado no lineal no ofrece una mejora significativa sobre el simple término lineal. Por lo tanto, el Modelo 2 es el preferido por ser más parsimonioso.
19 d) Análisis de la Salida del Modelo gam
La salida de summary.gam proporciona dos tablas de ANOVA: una para los términos paramétricos y otra para los no paramétricos.
summary(gam.m3)
Call: gam(formula = wage ~ s(year, 4) + s(age, 5) + education, data = Wage)
Deviance Residuals:
Min 1Q Median 3Q Max
-119.43 -19.70 -3.33 14.17 213.48
(Dispersion Parameter for gaussian family taken to be 1235.69)
Null Deviance: 5222086 on 2999 degrees of freedom
Residual Deviance: 3689770 on 2986 degrees of freedom
AIC: 29887.75
Number of Local Scoring Iterations: NA
Anova for Parametric Effects
Df Sum Sq Mean Sq F value Pr(>F)
s(year, 4) 1 27162 27162 21.981 2.877e-06 ***
s(age, 5) 1 195338 195338 158.081 < 2.2e-16 ***
education 4 1069726 267432 216.423 < 2.2e-16 ***
Residuals 2986 3689770 1236
---
Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
Anova for Nonparametric Effects
Npar Df Npar F Pr(F)
(Intercept)
s(year, 4) 3 1.086 0.3537
s(age, 5) 4 32.380 <2e-16 ***
education
---
Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1Interpretación del summary.gam(gam.m3): - Anova for Parametric Effects: Muestra la significancia de los términos lineales (education). El p-valor pequeño indica que education es un predictor significativo. - Anova for Nonparametric Effects: - La columna p-value(Linear) prueba la hipótesis \(H_0\): “La relación es lineal”. Para s(year, 4), el p-valor (0.348) es grande, sugiriendo que un efecto lineal es suficiente. - La columna p-value(Nonlinear) prueba la hipótesis \(H_0\): “La relación es nula”. Para s(age, 5), el p-valor (< 0.001) es pequeño, indicando que age tiene un efecto no lineal significativo.
20 e) Predicción con el Modelo GAM
Una vez seleccionado el mejor modelo (en este caso, gam.m2), se puede utilizar para realizar predicciones.
preds <- predict(gam.m2, newdata = Wage)
plot(preds, Wage$wage, xlab = "Valores Predichos", ylab = "Valores Observados", main="Predicciones vs. Observaciones")
abline(0, 1, col = "red", lwd = 2) # Línea de identidad (ajuste perfecto)
21 f) Uso de Regresión Local (LOESS) en GAM
En lugar de splines, se puede usar regresión local (lo()) como función de suavizado. El parámetro span controla el ancho de banda (el porcentaje de datos vecinos utilizados en cada ajuste local).
gam.lo <- gam(wage ~ s(year, df = 4) + lo(age, span = 0.7) + education, data = Wage)
# x11() # Descomentar para ver en una ventana separada
plot.Gam(gam.lo, se = TRUE, col = "green")
22 g) Introducción de un Término de Interacción No Paramétrico
Los GAM también pueden modelar interacciones entre predictores de forma no paramétrica, utilizando funciones de suavizado bivariadas (ej. lo(year, age)).
gam.lo.i <- gam(wage ~ s(year, df = 4) + lo(age, span = 0.7) + lo(year, age, span = 0.5) + education, data = Wage)Warning in general.wam(x, z, wz, fit$smooth, which, fit$smooth.frame, bf.maxit,
: general.wam convergence not obtained in 30 iterationssummary(gam.lo.i)
Call: gam(formula = wage ~ s(year, df = 4) + lo(age, span = 0.7) +
lo(year, age, span = 0.5) + education, data = Wage)
Deviance Residuals:
Min 1Q Median 3Q Max
-121.057 -19.810 -3.488 14.046 213.073
(Dispersion Parameter for gaussian family taken to be 1235.254)
Null Deviance: 5222086 on 2999 degrees of freedom
Residual Deviance: 3684803 on 2983.032 degrees of freedom
AIC: 29889.65
Number of Local Scoring Iterations: 3
Anova for Parametric Effects
Df Sum Sq Mean Sq F value Pr(>F)
s(year, df = 4) 1 22922 22922 18.557 1.702e-05 ***
lo(age, span = 0.7) 1 195729 195729 158.452 < 2.2e-16 ***
education 4 1074322 268581 217.429 < 2.2e-16 ***
Residuals 2983 3684803 1235
---
Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
Anova for Nonparametric Effects
Npar Df Npar F Pr(F)
(Intercept)
s(year, df = 4) 3.0 1.0418 0.37287
lo(age, span = 0.7) 1.2 4.9902 0.01958 *
lo(year, age, span = 0.5) 5.8 17.3006 < 2e-16 ***
education
---
Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1# La visualización de interacciones bivariadas produce un gráfico de contorno o perspectiva.
plot.Gam(gam.lo.i, se = TRUE, col = "orange")
Interpretación: La salida de summary y la prueba de ANOVA para los efectos no paramétricos indican que el término de interacción lo(year, age) es altamente significativo y no lineal, sugiriendo que la relación entre la edad y el salario cambia a lo largo de los años.